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Flux penetration into type-II superconducing films is simulated for transverse applied magnetic 
field and flux creep dynamics. The films contain macroscopic, non-conducting holes and we suggest a 
new method to introduce the holes in the simulation formalism. The method implies reconstruction 
of the magnetic field change inside the hole. We find that in the region between the hole and the 
edge the current density is compressed so that the flux density is slightly reduced, but the traffic 
of flux is significantly increased. The results are in good agreement with magneto-optical studies of 
flux distributions in YBa2Cu30^ films. 
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I. INTRODUCTION 

The behavior of vortex matter in superconductors can 
to a large degree be controlled by introducing artificial 
defects. It has been known for a long time that randomly 
distributed defects, created e.g. by neutron irradiation, 
allow a dramatic enhancement of the critical current den- 
sity, j c . One may reach more specific goals by tuning 
the arrangement of artificial defects. In particular, ex- 
periments on superconducting thin films have revealed 
a large number of interesting effects, including match- 
ing effects^ noise reduction in SQUIDs^ rectified vortex 
motion; 3 -^ anisotropy of and vortex guidance^ 

In parallel with the experimental progress, the theoret- 
ical understanding of how artificially created patterns in- 
teract with vortex matter is also developing. Interaction 
between a single vortex and a cylindrical cavity in a bulk 
superconductor was considered within the London ap- 
proximation in Ref. 0- This work extends the classical pa- 
per, Ref.[^, predicting the maximal number of flux quanta 
that can be trapped by a single hole. Current distribu- 
tion around a ID array of holes was calculated within the 
Ginsburg-Landau theory in Ref. 0. However, these theo- 
retical works consider a bulk superconductor, while most 
experiments are on patterned thin filmsJ*2£££& More- 
over, a realistic model should take into account the strong 
pinning of vortices in the superconducting areas around 
the artificial defects. When the defect size is much larger 
than the London penetration depth, one can consider the 
average vortex density B rather than individual vortices. 
Such an approach was used in Ref. [9j to simulate flux 
penetration into a thin film with a 2D array of holes. It 
allowed to explain an asymmetrical flux penetration due 
to asymmetry in the hole shape. At the same time, the 
case of an individual hole in a thin film has not yet been 
carefully analysed. A main purpose of the present work 
is to acquire details of flux and current distributions in a 
superconducting strip with one individual hole. 

An approximate picture of the current distribution 
around a non-conducting hole can be obtained within 
Bean's critical state models In the Bean model current 
stream lines are added from the edge with equal spac- 



ing representing the critical current density. The pres- 
ence of a hole forces the current to flow around it and 
hence pushes the flux front deeper into the sample. Both 
holes and sample corners give rise to so-called d-lincs 
where the current changes direction discontinuously.— 
They are seen as dark lines^ in images showing mag- 
netic flux distributions ! 13 ' 14 For example, 90° corners 
give 45° straight (i-lines^ while semicircular indentations 
of the edge give parabolic (i-lines J£ The magneto-optical 
image of Fig. [1] shows (i-lines spreading out from a cir- 
cular hole towards the flux-free region. The same hole 
also introduces another pattern: a darkened region start- 
ing from the hole and extending towards the edge. This 
pattern is similar to the one observed by Eisenmenger 
et al. , Ref. |T3. The pattern does not fit with the com- 
mon interpretation of the Bean model, which leaves the 
the currents between the hole and the edge unperturbed. 
Ref. [I?] discusses how to reinterpret the Bean model and 
explain the observed pattern as a second parabolic d-line. 
In this work, we will go further and do full dynamical sim- 
ulations of flux penetration taking into account the non- 
local electrodynamics of films as well as flux creep. Our 
results provide details of flux and current distributions 
in the vicinity of a hole and suggest a new interpretation 
for the observed anomaly. 



II. MODEL 

A. Single-connected superconductors 

Consider a type-II superconducting thin film placed 
in an increasing transverse magnetic field. The super- 
conductor responds by generating screening currents to 
shield its interior. The current density is highest at the 
edges where the Lorentz force eventually overcomes the 
pinning force, leading to penetration of flux. According 
to the Bean model, the vortices move only when the lo- 
cal current density exceeds the critical value, j c . A more 
realistic model for flux penetration also allows for flux 
creep at j < j c . Macroscopically, flux creep is introduced 
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FIG. 1: Left: Magneto-optical image of B z near a hole. Note 
the parabolic d-lines going upwards and a dark area going 
downwards from the hole. Right: a sketch of the strip with 
a circular hole indicating how peculiarities in the flux dis- 
tribution are related to bending of the current stream-lines. 
Notations: film half-width is a, distance from the edge to the 
hole center is s, and the hole radius is ro- 



through a highly non-linear current voltage relatio: 



E = pa 



Jj 



,15,18 



(1) 



where E is electric field, po a resistivity constant, j is 
current density, and n is the creep exponent. For thin 
films of YBa2Cu30 x , n is typically in the range from 10 
to 70 depending on temperature and pinning strength)! 9 . 

Flux dynamics of single-connected type-II supercon- 
ductors in transverse geometry has been described thor- 
oughly by E. H. Brandt. This work uses the same for- 
malism and hence we only give a short summary of the 
simulation basics, mainly following Refs. H5II20I and [2l|. 
The next section will be devoted to additional changes 
for multiply-connected samples. 

For films, it is a great simplification to work with the 

sheet current J(r) = J^^ 2 dz j(r, z), r = [x, y), in stead 
of the current density j. This is justified as long as thick- 
ness, d, is small compared to the in-plane dimensions but 
much larger than the London penetration depth, A. Fi- 
nite A can be handled with a small modification of the 
algorithm^ Since the current is conserved, V • J = 0, it 
can be expressed as J = V x zg, where g = g(r) is the 
local magnetization^ 

For single-connected thin films the Biot-Savart law can 
be formulated as 

B z (r, z)/no = H a + f d 2 r' Q(r, r', z) g(r'), (2) 

J A 

where H a is the applied magnetic field, and A is the 
sample area. The kernel Q represents the field generated 




FIG. 2: Simulated magnetic field distribution in a long strip, 
plotted in the style of magneto-optical images, where the in- 
tensity represents B z . At small applied field (left) the hole 
produces a field dipole and at large field (right) one can see 
the parabolic d-lines and a dark region between the hole and 
the edge, cf. the experimental image, Fig. [TJ ro/a — 0.1, 
s/a = 0.5, n = 19, H a /J c = 0.2 (left) and 1 (right), and 
PoH a = poJc/ad. 



by a dipole of unit strength^ 

1 2z 2 - (r - r') 2 



Q(r,r',2) = 



4-7T 



z 2 + (r — r') 



215/2 



(3) 



We discretize the kernel on an equidistant grid with grid 
points Ti and weights w and obtain 2 ^ 



}ij — Sij 



(4) 



where = l/47r|rj — | 3 for i ^ j and qu = 0. The 
function C depends on the sample geometry. It is given 
as 



C(r) 



dr r 



outside 



47r|r — r'l 3 



(5) 



The time evolution of g comes from the inverse of 
Eq. ©. 

g(r) = f d 2 r' Q' 1 ^^') [B z (r') - fx H a ], (6) 

J A 

where Q~ l for discrete problem is the matrix inverse of 
Eq. (0| . B z is given from Faraday's law as 

B z (r) = -(VxE), = V.(^-V S ), (7) 

with p = p \V g / J c \ n ~ 1 obtained from Eq. [TJ The right- 
hand side of Eq. is expressed only via g and H a so 
that time evolution of g can be found by integrating the 
equation numerically. 



3 



Bulk Bean model 



Film, n = 19 




FIG. 3: The current stream lines for the bulk Bean model 
(left) and for film with finite n (right). 



B. Superconductors with holes 

For macroscopic, arbitrarily shaped, single-connected, 
type-II superconducting films flux dynamics is fully de- 
scribed by Eq. ©. This basic equation can also be used 
for multiply connected samples, but in this case one needs 
to specify the dynamically changing value of g at the hole 
boundary. In Refs. [|| and this value was set to the 
lowest value of g along the hole perimeter. This method 
turned out to be quite feasible but unfortunately it can- 
not reproduce the discussed pattern of Fig. [1] Moreover, 
it also introduces unphysical net flux into the hole before 
the flux front has reached it. 

A completely different approach is to consider the holes 
as part of the sample, but ascribe to them a large Ohmic 
resistance or a strongly reduced J c 22- Then, Eq. © ap- 
plies to the whole sample including the holes, while the 
material law, Eqs. ([T]) and |(7J), is spatially non-uniform. 
This approach is physically justified but numerically chal- 
lenging due to huge electric field gradients. In addition, 
there still remain small but non-zero currents flowing 
within the holes. 

In this work we propose a new approach that does not 
require any additional assumptions, though requires a 
larger computational time. In this approach the integra- 
tion in Eq. ^ is extended over the whole sample area 
including the holes. Then the dynamics of g is described 
by the equation 

g(r) = f d 2 r> Q _1 (r,r') {B^^) + Bf'(r') - f , H a ], 

J A 

(8) 

where A is the sample area including the hole. Here we 



presented B z as a sum Bi ' 
>0) 



B^ 



where Bi' is nonzero 



only in the hole and B z is nonzero within the supercon- 
ducting areas. B z is calculated in the straightforward 
way using Eq. ||7J). The other term, B^\ is defined by 
two conditions. The first condition is that current does 
not flow beyond the superconducting areas, i.e., g is con- 
stant within the hole. This constant is determined by 
the second condition, that the total change of magnetic 
flux inside the hole is related to the electric field at its 



F. 




FIG. 4: Simulation results for a strip with a hole: the current 
stream lines (top), B z contour lines (middle), and E contour 
lines (bottom). Note that the electric field is greatly enhanced 
in the channel between the hole and the edgej^i H a / J c — 
0.3 and 0.9. The remaining parameters are the same as for 

pig. m 



boundary through Faraday's law, 



d 2 r B? 



hole 



tfl-E. 



(9) 



hole edge 



In order to find a B^ that satisfies the two conditions 
we use an iteration scheme. An initial guess, B z h ''°\ is 
substituted into Eq. flgj) to find g( h '°) inside the hole. The 
next approximation is found as 

Bj M) (r) = Bi h '°\r) - f d 2 r'Q(r, r')g < - h ' 0) (r') + K, 

J hole 

(10) 

where the constant K is chosen so that Eq. ^ is satis- 
fied. Bi h ' lS> is then inserted into Eq. (jSJ to find g( h,1 \ 
This gC 1,1 ) is in general non- uniform, but when the pro- 
cedure is repeated g^ h - n ^ becomes more uniform with ev- 
ery new iteration. A smart choice of the initial guess 
of £?i ?l ' '' is the final value at the previous time step, 



B 



M) 



B^ 



(r,t — At). With this choice only 



a couple of iterations are sufficient. 
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FIG. 5: Profiles of B z , J , g and E through y = for a strip 
with a hole. The curves correspond to applied fields H a / J c = 
0.1, 0.3, 0.5, 0.7, and 0.9, and E c = poJ c /d. The remaining 
parameters are the same as for Fig. [2] 



Note that the scheme presented here is in no way 
bound to the discrete formulation of the kernel, Eq. 
It can be used for any formulation as long as both the for- 
ward and inverse relations between g and B z are known. 
Further mathematical details are in appendix [X] 



III. STRIP WITH A CIRCULAR HOLE 

In this section, Eq. ([5]) is solved for an infinite su- 
perconducting strip in linearly increasing magnetic field. 
The strip is modeled using periodic boundary conditions 
in the y-direction, and examples of magnetic field distri- 
butions are given in Fig. O In the upper part one ob- 
serves regular flux penetration with maximum of B z at 




H a /J e 

FIG. 6: Total flux inside the hole, <& h = f. , d 2 r B,, as a 

' J hole ' 

function of H a , for various distances s from the edge. For low 
fields &h is zero since the flux front has not reached the hole 
yet. For high fields <E>h grows linearly with H a since the strip 
is saturated with J ~ J c . Hole radii are ro/a — 0.1. The 
remaining parameters are the same as for Fig. [2] 



the edges. Flux penetration in the lower half is strongly 
affected by the presence a small, circular, non-conducting 
hole. Note that the flux distribution is perturbed in a re- 
gion that significantly exceeds the hole dimensions. 

The left image of Fig. [2] corresponds to a small field 
for which the flux front has not reached the hole yet. In 
this case the hole shows up as a field dipole, i n ag reement 
with magneto-optical observations; cf. fiefs. [171 and l25l. 
Namely, there is positive field at the farther side of the 
hole and a negative field at the side closer to the film 
edge. The negative fields shrinks when the flux front 
reaches the hole, but the asymmetry of the flux distribu- 
tion inside the hole remains, as seen in the right image. 
As expected, the front becomes distorted so that the pen- 
etration is significantly deeper in the vicinity of the hole. 
For the full penetration image of Fig. [2] one also clearly 
see the e?-lines as dark line originating at the hole and 
directed towards the middle of the strip. Such d-lincs 
were first described in Ref. |Tl| within the Bean model 
framework and they arc called d-lines because current 
changes direction discontinuously there. The disconti- 
nuity is most clearly seen in current stream line plot of 
Fig. [3J (left). For the Bean model, d- lines from circular 
holes are parabolic and by convention d-lines from small 
holes inside superconductors are often called parabolas. 
In the presence of flux creep the change of current direc- 
tion is smeared as follows from Fig. [3J (right). However, 
the d-lines are still clearly visible, at least for n ^> 1. 

Comparing the two panels of Fig. [3J we notice a qual- 
itative difference between the current flow in the bulk 
Bean model and for films under the creep. In the Bean 
model the current density is everywhere constant and all 
the current that is blocked by the hole turns towards the 
strip center. The region between the hole and the edge 
is hence unaffected by the presence of the hole. For film 
creep dynamics this is no longer true and a certain frac- 
tion of the current will force its way here. As a result, 
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the current density is enhanced which is seen as denser 
stream lines in Fig. [3] (right). Since the stream lines bend 
they create the feature visible in the flux distribution of 
Fig. [5] a slightly darkened region starting at the hole 
and widening towards the edge. This feature can also be 
observed experimentally; cf. Fig. [TJ It was analysed in 
detail in Ref. [13 and interpreted in terms of the Bean 
model as additional parabolic d-lincs. Our experiment 
and simulations suggest a different interpretation. We 
believe that one should speak about an area of reduced 
flux density rather than new d-lines. Moreover, the ap- 
pearance of this area is due to locally enhanced current 
density, hence it cannot be explained within the Bean 
model, postulating J = J c . An enhanced current den- 
sity also implies a strongly enhanced electric field. This 
is clearly seen in Fig. 0] showing the contour lines of E. 
A locally enhanced E means that there is an exception- 
ally intensive traffic of magnetic flux through the channel 
between the edge and the hole. The channel width is ap- 
proximately given by the hole diameter, but increases 
slightly towards the edge. The width depends in gen- 
eral on the distance to the edge and the creep exponent 
n. Both larger distance and smaller n tend to make the 
channel wider. 

After arrival to the hole, the flux is further directed in 
the fan-shaped region between the d-lines. Electric field 
within this region is also relatively high, again implying 
an intensive flux traffic. This situation is similar to the 
case of a semicircular indentation at the sample edge con- 
sidered in Refs. [r| [H, and H3- The hole thus strongly 
rearranges trajectories of flux flow. 

The above discussion is further confirmed by profiles 
of B z , J, g, and E through the line y = shown in 
Fig. The J profiles show features commonly observed 
in strips^ i.e., plateaus with values ~ J c in the pene- 
trated regions and shielding currents with J < J c in the 
Mcissncr regions. The profiles show clearly the enhanced 
J and E between the edge and the hole. It is also in- 
teresting to see the negative B z for low values and how 
the negative values gradually vanish when the main flux 
front gets in contact. 

Fig. [5] shows the total flux in a circular hole, = 
J hole d 2 r B z , as a function of the applied field H a for 
various distances between the hole and the edge. In the 
beginning, <&h w 0, until the main flux front is in contact 
with the hole. Then it starts to increase. For high fields 
<5>h grows almost linearly with H a at a universal growth 
rate determined by the hole area. The linear rate is not 
just the case for small holes in strips, but has also been 
found for e.g. ring geometry^ Note that for small fields 
<f>/j is close to, but not exactly zero. The reason is the cre- 
ation of two additional flux fronts: one positive towards 
the flux-free region and one negative towards the edge, 
as also seen in Fig. [5] Only when integrating B z over a 
larger area that includes this additional penetration one 
finds that the total flux is exactly zero^ This integral 
is also a good consistency check of the boundary condi- 
tion implementation, since a wrong value of g^ h ' tends to 





FIG. 7: A square with two slits (only the lower half is shown). 
Top: Sample sketch, experimental magneto-optical image of 
YBa2Cu30;c film, and simulated magnetic field distribution. 
Bottom: current stream lines, B z and E contour lines at 
H a /Jc — 0.3 and 0.9, with n = 19 and fJ,oH a = poJc/ad. 
Note the strongly enhanced J and E between the slit and the 
edge and the complicated set of d-lines at full penetration. 

introduce a net, unphysical flux in the hole. 

IV. SQUARE WITH TWO SLITS 

This section presents results of simulations of a square 
superconducting film with rectangular slits. The sam- 
ple geometry is chosen to reassemble one particular 
YBa2Cu30 a . film and comparison of the simulation with 
a magneto-optical image of the sample is shown in Fig. [7] 
The experimental film thickness is 250 nm, and side 
lengths are 2.5 mm. The two slits have been cut out 
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with a laser. Details of the film preparation can be found 
elsewhere 4i 

The experiment and simulations show a great similar- 
ity both in large and in the details. The flux density is 
considerably enhanced everywhere along the slit edges, 
and reaches the maximal values at the upper corners. 
Our main result found for circular holes holds true also 
for rectangular slits. Namely, we again find a distinct 
dark region starting at the slit and widening towards the 
edge. It can be attributed to the over-critical current 
density in that region, which is clearly seen in the cur- 
rent stream-line plot. A new result for slits is a slightly 
brightened regions near the upper corners that appear 
due to concave current turns. A similar situation arises 
in superconductors of some other shapes having concave 
corners, e.g. in crosses^ 

There also exist a few minor discrepancies between flux 
distributions obtained in the simulations and in the ex- 
periment of Fig. [7J The most notable is the details of the 
region of reduced B z at the side of slits close to the edge. 
The values of B z appear to be less in the simulation than 
in experiment. This might be caused by simplifications, 
like the disregarded i?-dcpendcncy of the material law or 
the simplification of using the sheet current in stead of 
the true current density. 

V. SUMMARY 

We have proposed a new method for treating boundary 
conditions of non-conducting holes inside macroscopic, 
typc-II superconducting films. The key point is to recon- 
struct the at first unknown B z inside the holes, at each 
time step of the simulation. The method is capable of 
handling any number of holes of arbitrary shape. 

The simulations of flux dynamics assuming a material 
law E ~ j n reproduce very well flux distributions ob- 
served by magneto-optical imaging in YB^CusOa; films, 
for circular holes as well as rectangular slits. In par- 
ticular, they demonstrate a significant enhancement of 
current density in the region between a hole (slit) and 
the edge leading to a more intensive traffic of flux. This 
region appears darker in magneto-optical images due to 
a slight bending of current stream lines. 

We thank C. Romcro-Salazar and Ch. Jooss for fruit- 
ful discussions and M. Baziljevich for experimental data 
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Norwegian Research Council, Grant No. 158518/431 
(NANOMAT) and by FUNMAT@UIO. 



APPENDIX A: NUMERICAL DETAILS 

The simulations are run on a TV x N square grid. 
The creep exponent and the ramp rate are n = 19 and 
jj,oH a = poJ c /ad, a regime in which creep is low, but 
not negligible. Changing n would only do quantitative 
changes to the results. For small exponents the plateaus 
of current profiles, like Fig. [5J would be less flat and there 
would also be more current compressed between the holes 
and the edge. 

The main limiting factor of the simulations is mem- 
ory consumption since the kernel matrix Q, Eq. ([5]), 
has dimension iV 2 x TV 2 . The simulations are run with 
TV = 100 grid points, which yields a kernel matrix of 
dimension 5000 x 5000, when the sample symmetry has 
been exploited^ 

The kernel Q in Eq. ^ depends explicitly on the sam- 
ple shape. Since the strip is infinite in the y-direction, 
Q should be computed via an infinite sum over strip seg- 
ments. However a good approximation is achieved with 
only one segment on each side of the "main" strip. The 
strip segments further away contain zero net current and 
the dipole like character means that they have a neg- 
ligible effect. A good accuracy of this approximation 
was checked by comparing the Meissner state width, 6, 
obtained for very high n with the analytical film Bean- 
model result^ b = a/ cosh(irH a / J c ). 

The reconstruction of B z inside the hole, Eq. (|10p . 
need not use the full Q from Eq. The best is to 

use a smaller kernel, Q, also generated with Eq. ([4]), but 
only including points inside the hole. Fast convergence 
of Eq. (JTUJ) is achieve by ignoring currents at the hole 
perimeter, which means that Q should use C(r) = 0. 

The most difficult numerical problem in our method 
is the calculation of the electric field at the boundary 
in Eq. @. The electric field is given by the power law, 
Eq. ([T]), and is largely fluctuating between neighboring 
grid points. A stable way to handle this is to take the 
average of only the most significant values of E and use 
2nr and irr 2 for the hole circumference and area. 
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